Force indeterminacy in the jammed state of hard disks 



Tamas Unger 

Dept. of Theoretical Physics, Budapest University of Technology and Economics, H-llll Budapest, Hungary and 
Institute of Physics, University Duisburg-Essen, D-47048, Duisburg, Germany 

Janos Kertesz 

Dept. of Theoretical Physics, Budapest University of Technology and Economics, H-llll Budapest, Hungary 

Dietrich E. Wolf 

Institute of Physics, University Duisburg-Essen, D-47048, Duisburg, Germany 
(Dated: February 2, 2008) 

Granular packings of hard discs are investigated by means of contact dynamics which is an appro- 
priate technique to explore the allowed force-realizations in the space of contact forces. Configura- 
tions are generated for given values of the friction coefficient, and then an ensemble of equilibrium 
forces is found for fixed contacts. We study the force fluctuations within this ensemble. In the 
limit of zero friction the fluctuations vanish in accordance with the isostaticity of the packing. The 
magnitude of the fluctuations has a non-monotonous friction dependence. The increase for small 
friction can be attributed to the opening of the angle of the Coulomb cone, while the decrease as 
friction increases is due to the reduction of connectivity of the contact-network, leading to local, 
independent clusters of indeterminacy. We discuss the relevance of indeterminacy to packings of 
deformable particles and to the mechanical response properties. 



Jamming [l| has been in the focus of recent studies 
because it occurs in a great variety of phenomena like 
structural and spin glasses, colloidal systems, vehicular 
traffic and granular media. The characterization of the 
jammed state is therefore crucial and can perhaps be best 
achieved in granular systems. Many intriguing proper- 
ties of granular packings originate from the microscopic 
force transmission through a contact structure, where 
non-linearity and disorder are known to be crucial. It 
is an essential but not resolved question how the highly 
inhomogeneous force-network influences the macroscopic 
stress transmission in dense granular media. 

Since the deformations of the grains are usually much 
smaller than their size, a very useful reference system for 
granular matter is that of rigid (undcformable) particles 
[H y. Ll [j . It is known that random packings of frictional 
rigid disks or spheres exhibit a hyperstatic structure H 
la, 0] : the number of the linear equilibrium equations of 
the grains, which relate the unknown contact forces to 
the external load, is too small to determine the contact 
forces uniquely. Therefore many mechanically admissible 
force-networks are possible in the same packing geometry 
and for the same external load, which define an ensemble 
of force-configurations. 

This ensemble recentl y h as received much attention 
H H H El El El El IH due to the idea that some 
macroscopic properties of jammed granular systems can 
be derived based on an ensemble average over the admis- 
sible force-states The determination of force distri- 
bution in or Green function in E3 are based on this 
approach. 

Another interesting aspect of the force-ensemble is re- 
lated to the behavior of the system under external per- 
turbations. Packing structures where contact forces are 
unique or strongly restricted appear to be fragile: slight 



change of the load can cause rearrangements of the parti- 
cles [3, 0| . The question arises whether a packing that 
exhibits many possible realizations of equilibrium forces 
becomes more robust against perturbations. 

The results of this Letter provide nontrivial informa- 
tion also for packings of deformable particles: The actual 
network of contact forces (which is uniquely determined 
by the elastic deformations) must be contained in the 
force-ensemble calculated for the same contact geometry 
assuming the particles (in their deformed shape) would 
be perfectly rigid. Moreover, for a finite system of suffi- 
ciently rigid particles the contact geometry can be arbi- 
trarily close to the ideal one obtained for perfect rigidity. 
Which of the solutions in the force ensemble is realized, 
depends e.g. on the elasticity of the individual grains. 
Here we address the question, how strong the restrictions 
provided by the force ensemble are. 

Again another but closely related issue is that of hard 
particle simulations, where the dynamics is seemingly 
ambiguous due to the indeterminacy of forces |l3j . 

The above problems indicate the significance of the 
force-ensemble, however very little is known about its 
properties. In this Letter some characteristics of the en- 
semble are revealed, where emphasis is put on the influ- 
ence of friction. 

In the recent literature H it was suggested that 
all elements in the ensemble of admissible force config- 
urations are realized with equal probability. This mi- 
crocanonical approach can be regarded as a restricted 
version of the Edwards ensemble H 0. In the 
following we also address the validity of this assumption. 

Let us consider n rigid, cohesionless disks. A configu- 
ration of the contact forces {Fj} (where i is the contact 
index) is called admissible or a solution if two conditions 
are fulfilled: the equilibrium and the Coulomb conditions. 
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The first one requires force and torque balance at each 
grain, while the Coulomb condition reads: 

|(Fj) t |</i(Fj)„ (1) 

for the normal and tangential force at each contact, where 
fi is the friction coefficient. For /i > no additional 
condition is needed to exclude tensile forces. 

Next we show that the solutions form a convex set. 
The space of contact forces T is defined (for fixed con- 
tact network) as an N c x d dimensional vector space, 
where each point represents a force-configuration {F;}. 
iV c is the number of contacts, and d the space dimen- 
sion (i.e. each contact force component represents one 
degree of freedom). Let S be the subset of admissible 
states in T under some fixed external forces. For a reg- 
ular packing of disks S is known to be a convex polyhe- 
dron [TH but it is easy to see that convexity is satisfied 
in any case: shape of the particles, disorder, dimension- 
ality or friction do not matter. Convexity means that 
if {Fi} and {Fj + AFJ are solutions then {F. t + AAFJ 
is a solution as well for < A < 1. First, the equi- 
librium condition holds: Both given force-configurations 
provide equilibrium against the external load, thus their 
difference {AFj} corresponds to zero load and exerts no 
total force or torque on the particles. Therefore it can be 
scaled freely (unrestricted A) and added to an admissible 
state, that does not violate the linear equilibrium equa- 
tions. Second, the Coulomb condition is satisfied simply 
because for each contact i the d-dimcnsional Coulomb 
"cone" is a convex set and therefore must contain the 
component Fj + A AFj , with < A < 1 . 

The solution set S reflects basically the properties of 
the contact-network, therefore when studying S it is cru- 
cial what kind of packing structure is considered. In 
real processes which lead to jamming, the microscopic 
structure is not prescribed but develops spontaneously 
up to the point, where further rearrangements against 
outer driving forces are blocked. This self-organized tex- 
ture is an important feature of granular materials [l5j 
which is disregarded in models using, e.g., regular ar- 
rangements 0, • Therefore the packings studied be- 
low were constructed with discrete element simulations 
where the particles obeying Newton's dynamics build up 
the contact-network in a compression process. In these 
jammed configurations we search for various solutions of 
the contact forces and study the influence of friction on 
the properties of S. 

A detailed description of our method of con- 
structing the packings and exploring admissible force- 
configurations can be found in |6|], here only a short 
review is given. With the help of the contact dynam- 
ics algorithm [H E3 a 2D system of 200 rigid disks is 
compressed along the vertical axis between two horizon- 
tal plates. Horizontally periodic boundary conditions are 
applied, gravity is set to zero, disk radii are uniformly dis- 
tributed between R and 2R, the horizontal system width 
is 42i?. We wait till the packing jams (relaxes into equi- 
librium) under the constant force of compression. Then, 



to avoid the effect of the straight plates, only the mid- 
dle part of the static configuration is considered for fur- 
ther investigation: this is a horizontal slice of height 28R 
throughout the whole width in the bulk away from the 
plates. We retain the contact forces at the top and bot- 
tom perimeter of the slice as fixed boundary forces, thus 
they provide the external load on the system. The plates 
and the disks outside the slice can be left away. 

After that the exploration of the admissible force- 
solutions follows for this fixed arrangement of disks. We 
start with the force state that appeared at the jam- 
ming and perturb all contact forces randomly |24| , which 
leads out of equilibrium and violates the Coulomb con- 
dition. This perturbed state serves as the input for the 
Gauss-Seidel-like iterative solver of the contact dynam- 
ics method. This iterative algorithm lets the forces relax 
into a consistent state, providing a (possibly) new so- 
lution 0, Il8| . The perturbation and relaxation can be 
repeated many times always starting from the last solu- 
tion (a kind of random walk in the force space); in that 
way it is possible to sample points from S. 

Based on this collection of force solutions we can assess 
the differences between admissible states and study the 
problem of force indeterminacy. The main feature of S 
that we found in these self-organized structures is that 
the admissible force-networks are rather similar: The 
pattern of strong force lines changes little from one re- 
alization to the other, showing that the contact-network 
imposes strong restrictions on the force-configuration. 

For each contact force Fj its variance (SFi) 2 is calcu- 
lated over the measured realizations. The ratio 

r? = <^>/<|F|) (2) 

represents the ensemble fluctuation in 5, thus it can be 
regarded as a measure of ambiguity of the forces. (•} 
means the average over all contacts. The force ambiguity 
n has to be distinguished from the degree of indetermi- 
nacy which refers to the dimension of the affine subspace 
of force configurations solving the equilibrium conditions 
(without the restrictions due to the Coulomb cones). 

To investigate the effect of friction a new packing 
is constructed for each value of [i before sampling the 
solutions. The force ambiguity r\ is plotted in Fig. ^ 
(full circles). Values of n around reflect the ac- 

curacy level of our calculation and the corresponding 
force-configurations can be regarded as identical with 
this tolerance. In the zero friction limit the force am- 
biguity disappears confirming isostaticity of frictionless 
packings 0, 0, H3- For small \i the force ambiguity 
grows proportionally with friction, however for larger fi 
it decreases again. The largest ambiguity of the forces 
is found around fi « 0.1. Despite the further opening of 
the Coulomb angle fluctuations are getting smaller, even 
fully determined states are found for strong friction. 

The behavior of r\ results from two competing effects: 
first, increasing friction provides larger freedom locally 
for the tangential forces, second, it also stabilizes the 
system in a less dense state [2l| causing lower connectiv- 
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FIG. 1: Force ambiguity 77 (full circles) and average coordi- 
nation number 2 (open circles) as functions of the friction 
coefficient p. For comparison, squares connected by line show 
the 77 values for a configuration of disks that was constructed 
without friction. 

ity of the contact-network (open circles in Fig.^l, which 
reduces force ambiguity. One can separate the two effects 
by fixing the configuration and letting the Coulomb an- 
gle alone influence 77: We generated one packing without 
friction but switched on friction before sampling force- 
configurations. The results obtained this way (squares in 
Fig. ^) provide monotonously increasing fluctuations, as 
expected. Compared to the original data (full circles) de- 
viations appear only on the right side of the figure, where 
the changes in the connectivity become important, while 
the behavior on the left side is governed by the first ef- 
fect. For small p the average coordination number of the 
configuration is essentially the same as in the friction- 
less case, where from isostaticity N c w 2n follows. This 
gives us the degree of indeterminacy: 2N C — 3n « N c /2, 
since there are two unknown force components per con- 
tact and three equations per disk due to force and torque 
balance. Thus we conclude that for tiny friction there is 
a small but high- dimensional set of force-solutions in the 
2N C dimensional force-space, and its size goes to zero 
with vanishing friction. Similarly for spheres in three di- 
mensions one obtains an ./V c -dimensional solution set S 
within a 3A^ c -dimensional force space T . 

For large p the dimension of S is strongly reduced due 
to the decreasing number of contacts. In our small system 
we found that dim(iS) can reach even zero, allowing only 
one single force- configuration. This case corresponds to 
the marginal rigidity state found in experiments . 

The regression of the degrees of freedom occurs in 
an interesting way: the indeterminacy gets localized in 
space into small subgraphs of the contact-network, which 
are surrounded by determined forces, i.e. a relatively 
large ambiguity is present but only in a small part of 
the system (Fig. 0b). The pattern of the fluctuation- 
bearing contacts can be visualized by plotting the dif- 
ference between any two admissible force-configurations. 
We found the same subgraphs as in Fig. 0b also for 
other arrangements of boundary forces, showing that this 




FIG. 2: (Color online) The difference between two admissible 
force- networks for (a) p = 0.1 (b) p = 0.5. Only normal force 
differences are indicated with different colors depending on 
their sign. 

indeterminacy-pattern is indeed a property of the pack- 
ing texture. Each of the two subgraphs shown in Fig. 0b 
is statically indeterminate, carries only one degree of free- 
dom and cannot be reduced further because the deletion 
of one particle or one contact would cancel the inter- 
nal indeterminacy. We call such subgraphs elementary 
clusters. They can be regarded as geometric units of in- 
determinacy. 

If the connectivity is high the formation of elemen- 
tary clusters is more probable, which suggest the follow- 
ing picture: For small friction many overlapping elemen- 
tary clusters are formed so that two admissible solutions 
generically differ throughout the system (Fig. 0a). As 
N c is reduced the density of the elementary clusters p 
decreases and the indeterminacy gets localized into small 
separated domains. Around p — 1 the density p becomes 
so small that only a few elementary clusters are present 
due to the finite system size. This explains the strong 
scattering of the data for 77 in Fig. 

The spatial localization raises the question of a per- 
colation transition. In case of small p the separated 
domains carry force-fluctuations independently of each 
other, therefore we think that 77 becomes a well defined 
intensive quantity for large systems. However if the inde- 
terminacy percolates through the system the overlapping 
elementary clusters provide fluctuating boundary forces 
for each other, thus the indeterminacy of forces is en- 
hanced with growing system size. Simulations up to 500 
particles show this size dependence, but it is not clear 
what happens in the thermodynamic limit. 

Finally we investigate the dynamically created force- 
configuration {F^o}, which is determined by the con- 
struction history. Our findings indicate that this state 
is more "central" than typical points in the solution set: 
We generate 20 initial configurations with p = 0.01 and 
sample for each of them 100 points randomly in S. Their 
(vectorial) average is regarded as the center of S. Then 
we measure the Euclidean distances I of the sampled 
points from the center. The histogram of the distances in 
units of their average I is shown in Fig. together with 
the histogram of the distances £q of the initial, dynam- 
ically generated 20 points from the centers of the cor- 
responding sets S. The two histograms clearly indicate 
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FIG. 3: Histograms of the distributions of normalized dis- 
tances of dynamically generated (full circles) and randomly 
sampled (open circles) points in the sets S for fi = 0.01. The 
inset shows a two dimensional cross section of a high dimen- 
sional solution set. The dynamically constructed force state is 
marked by the arrow. The white area belongs to S, while out- 
side S the gray scale indicates the violation of the Coulomb 
condition (darker means smaller violation). 

that the initial points are closer to the center on aver- 
age than the randomly sampled ones. Assuming that the 
distribution of the random sampling of S is close to a uni- 
form one, we conclude that the force configurations of the 
dynamically generated jammed states are not uniformly 
distributed in the set S. 

That the original force configuration is "closer to the 
center" is not in contradiction to the fact that we always 
find it at the edge of two dimensional cross sections of the 
high dimensional solution set S (see inset of Fig. [SI. We 



suggest the following physical picture: A contact with 
large mobilization of friction (Ft/fiF n « 1) is less stable 
against perturbations. Near the end of the relaxation 
process small collisions "shake" the already established 
contacts reducing the possibility that the contact remains 
on the verge of sliding. However, the system comes to 
rest finally by the marginal fulfillment of the Coulomb 
criterion at some contacts. 

Our results show a significant difference between dis- 
tributions of the solutions sampled by the random walks 
plus relaxation and of those relaxed physically. The uni- 
formity of the (unbiased) random walk based sampling 
cannot be proved due to the high dimensionality of the 
problem, however, the distance distribution of the points 
should be rather robust just because of this high dimen- 
sionality. Therefore we consider the observed discrepancy 
though not as a proof but as a strong indication of the 
violation of the microcanonical assumption for the phys- 
ically realized solutions. 

It is expected that the ambiguity of forces for a given 
geometry has implications for the mechanical behavior. 
We regard the following preliminary result as an indica- 
tion of such an effect. For a horizontal layer of hard disks 
settled under gravity we applied a point-force downwards 
on the free surface, just strong enough to cause local rear- 
rangement. We measured the depth of the rearrangement 
zone and obtained non-monotonous dependence on /i: It 
is larger for small and large friction coefficients, and has a 
minimum at /i s=s 0.1, right where n reaches its maximum. 
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